#define NULL 0x0

#define BITS_IN_INT_GEN 32

#if (BITS_IN_INT_GEN == 32)
#define INT_MOD_MASK 0xffffffff
#define FLT_MULT (0.25/(unsigned int)(1<<30))
#else
#define INT_MOD_MASK ((unsigned)(1<<BITS_IN_INT_GEN)-1)
#define FLT_MULT (1.0/(1<<BITS_IN_INT_GEN))
#endif

#define INT_MASK ((unsigned int)INT_MOD_MASK>>1)
#define MAX_BIT_INT (BITS_IN_INT_GEN-2)
#define INTX2_MASK ((1<<MAX_BIT_INT)-1)

#define RUNUP (2*BITS_IN_INT_GEN)

#define GS0 0x372f05ac

#define BASE_ADDR_NSTART(thread_num, gennum, genptr) (0 + (gennum) + (__global unsigned int*)genptr) 
#define BASE_ADDR_R0(thread_num, gennum, genptr) ((thread_num) * 1280  +  (__global unsigned int*)(genptr))
#define BASE_ADDR_R1(thread_num, gennum, genptr) ((thread_num) * 1280 * 2 + (__global unsigned int*)(genptr))
#define BASE_ADDR_SI(thread_num, gennum, genptr) ((thread_num) * 1280 * 3 + (gennum) + (__global unsigned int*)(genptr))
#define BASE_ADDR_HPTR(thread_num, gennum, genptr) ((thread_num) * 1280 * 4 + (gennum) + (__global int*)(genptr))

#define MAP(I, thread_num) ((I) * (thread_num))

#define MAX_THREAD_PERBLK_NUM 512  //TEMP AND NOT REASONABLE, TO AVOID

//#define DEBUG

#define HPTR_SIZE 64

/*************************************************************************/
/*************************************************************************/
/*                  STRUCTURES AND GLOBAL DATA                           */
/*************************************************************************/
/*************************************************************************/

typedef enum _CLSprngErr
{
  kCLSprngErrNoError = 0,
  kCLSprngErrGenNumExceedsTotal,
  kCLSprngErrTotalGenIsNegative,
  kCLSprngErrGenNumIsNegative,
  kCLSprngErrUnknownParameter
}CLSprngErr;

typedef struct _rngen 
{
  int globalSeed;
  int numGenerators;
  unsigned int si[1280];
  unsigned int r0[1280];
  unsigned int r1[1280];
  unsigned int nstart[1280];
  int seed;
  int init_seed;
  unsigned int lValue;
  unsigned int kValue;
  unsigned int LSBS;
  unsigned int first;
  int total_gen;
  int hptr;
}rngen;

//Move constants local->global
#if 0
void clSPRNG_packConstants(__global rngen* generator,
                           __local rngen_const* constants)
{
  barrier(CLK_LOCAL_MEM_FENCE);

  if(get_local_id(0) == 0)
  {

    generator->constants.seed = constants->seed;
    generator->constants.init_seed = constants->init_seed;
    generator->constants.lValue = constants->lValue;
    generator->constants.kValue = constants->kValue;
    generator->constants.LSBS = constants->LSBS;
    generator->constants.first = constants->first;
    generator->constants.total_gen = constants->total_gen;

    for(unsigned int i = 0; i < HPTR_SIZE; i++)
    {
      generator->constants.hptr[i] = constants->hptr[i];
    }
  }

  barrier(CLK_LOCAL_MEM_FENCE);
}

//Move contants from global->local
void clSPRNG_unpackConstants(__global rngen* generator,
                             __local rngen_const* constants)
{
  barrier(CLK_LOCAL_MEM_FENCE);

  if(get_local_id(0) == 0)
  {
    constants->seed = generator->constants.seed;
    constants->init_seed = generator->constants.init_seed;
    constants->lValue = generator->constants.lValue;
    constants->kValue = generator->constants.kValue;
    constants->LSBS = generator->constants.LSBS;
    constants->first = generator->constants.first;
    constants->total_gen = generator->constants.total_gen;

    for(unsigned int i = 0; i < HPTR_SIZE; i++)
    {
      constants->hptr[i] = generator->constants.hptr[i];
    }
  }

  barrier(CLK_LOCAL_MEM_FENCE);
}
#endif

#define kNumParams 11

__constant unsigned int validL[kNumParams] = 
  {1279, 17, 31, 55, 64, 127, 521, 521, 607, 607, 1279};  

__constant unsigned int validK[kNumParams] = 
  {861, 5, 6, 24, 31, 97, 353, 168, 334, 273, 418};

__constant unsigned int validLSBS[kNumParams] =
  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};

__constant unsigned int validFirst[kNumParams] =
  {233, 10, 2, 11, 14, 21, 100, 83, 166, 105, 208};

__kernel void clSPRNG_getGeneratorSize(__global unsigned int* size)
{
  if(get_global_id(0) == 0)
  {
    *size = sizeof(rngen);
  }
}


/*************************************************************************/
/*************************************************************************/
/*                    ERROR PRINTING FUNCTION                            */
/*************************************************************************/
/*************************************************************************/

void errprint(__constant char *level, 
              __constant char *routine, 
              __constant char *error)
{
}

/*************************************************************************/
/*************************************************************************/
/*            ROUTINES USED TO CREATE GENERATOR FILLS                    */
/*************************************************************************/
/*************************************************************************/

/**************************/
/* function bitcnt:       */
/**************************/
int bitcnt(int x)
{
  unsigned int i = 0, y;

  for (y=(unsigned int)x; y; y &= (y-1))
  { 
    i++;
  }

  return(i);
}

/**************************/
/* function advance_reg:  */
/**************************/
void advance_reg(int *reg_fill)
{
  const unsigned int mask = 0x1b;
  __local int adv_64[4][2];

  int i,new_fill[2];
  unsigned int temp;

  adv_64[0][0] = 0xb0000000;
  adv_64[0][1] = 0x1b;
  adv_64[1][0] = 0x60000000;
  adv_64[1][1] = 0x2d;
  adv_64[2][0] = 0xc0000000;
  adv_64[2][1] = 0x5a;
  adv_64[3][0] = 0x80000000;
  adv_64[3][1] = 0xaf;
  new_fill[1] = new_fill[0] = 0;
  temp = mask<<27;

  for (i = 27; i >= 0; i--) 
  {
    new_fill[0] = (new_fill[0] << 1) | (1 & bitcnt(reg_fill[0] & temp));
    new_fill[1] = (new_fill[1] << 1) | (1 & bitcnt(reg_fill[1] & temp));
    temp >>= 1;
  }

  for (i = 28; i < 32; i++) 
  {
    temp = bitcnt(reg_fill[0] & (mask << i));
    temp ^= bitcnt(reg_fill[1] & (mask >> (32 - i)));
    new_fill[0] |= (1 & temp) << i;
    temp = bitcnt(reg_fill[0] & adv_64[i-28][0]);
    temp ^= bitcnt(reg_fill[1] & adv_64[i-28][1]);
    new_fill[1] |= (1 & temp) << i;
  }

  reg_fill[0] = new_fill[0];
  reg_fill[1] = new_fill[1];
}

/**************************/
/*   function get_fill:   */
/**************************/

int get_fill(__global unsigned int* n, 
             __global unsigned int* r, 
             __global rngen* generator)
{
  int i;
  int j;
  int k;
  int temp[2];
  int length;
  int total_gen = generator->total_gen;

  length = generator->lValue;

  temp[1] = temp[0] = n[0] ^ generator->seed;

  if (!temp[0])
  {
    temp[0] = GS0;
  }

  advance_reg(temp);
  advance_reg(temp);

  r[0] = (INT_MASK & n[0]) << 1;

  for (i = 1; i < length - 1; i++) 
  {
    advance_reg(temp);
    r[i] = (INT_MASK & (temp[0] ^ n[i])) << 1;
  }

  r[length - 1] = 0;

  k = generator->first + generator->LSBS;

  for (j = generator->first; j < k; j++)
  {
    r[j] |= 1;
  }

  return 0;
}

/*************************************************************************/
/*************************************************************************/
/*            SI_DOUBLE: updates index for next spawning                 */
/*************************************************************************/
/*************************************************************************/

void si_double(__global unsigned *a,  
               __global unsigned *b, 
               int length, 
               int total_gen)
{
  int i;

  if (b[MAP(length - 2, total_gen)] & (1 << MAX_BIT_INT))
  {
    //error;
  }
  a[MAP(length - 2, total_gen)] = (INTX2_MASK & b[MAP(length - 2, total_gen)]) << 1;

  for (i = length - 3; i >= 0; i--) 
  {
    if (b[MAP(i, total_gen)] & (1 << MAX_BIT_INT)) 
      a[MAP(i + 1, total_gen)]++;
    a[MAP(i, total_gen)] = (INTX2_MASK & b[MAP(i, total_gen)]) << 1;
  }
}

/*************************************************************************/
/*************************************************************************/
/*            GET_RN: returns generated random number                    */
/*************************************************************************/
/*************************************************************************/

int get_rn_int(__global rngen *generator, 
               int flag)
{
  unsigned int new_val;
  int lptr;
  int hptr = generator->hptr;
  int lval, kval;
  int total_gen = generator->total_gen;
  unsigned int tmp, tmp0, tmp1;
  unsigned int lValue = generator->lValue;
  unsigned int kValue = generator->kValue;

  lptr = (hptr + kValue) % lValue;

  if(flag != 0)
  {
    tmp0 = INT_MOD_MASK & 
      (generator->r0[hptr] + generator->r0[lptr]);
    tmp1 = INT_MOD_MASK & 
      (generator->r1[hptr] + generator->r1[lptr]);

    new_val = (tmp1 & (~1)) ^ (tmp0 >> 1);
    
    generator->r0[hptr] = tmp0;
    generator->r1[hptr] = tmp1;

    hptr = --hptr < 0 ? lValue - 1 : hptr;
    hptr = --hptr < 0 ? lValue - 1 : hptr;
    
    generator->hptr = hptr;
  }
  else
  {
    tmp = (lptr + kValue) % lValue;

    tmp0 = INT_MOD_MASK & 
      (generator->r0[tmp] + generator->r0[lptr]);
    tmp1 = INT_MOD_MASK & 
      (generator->r1[tmp] + generator->r1[lptr]);

    generator->r0[lptr] = tmp0;
    generator->r1[lptr] = tmp1;

    tmp0 = INT_MOD_MASK & (generator->r0[hptr] + tmp0);
    tmp1 = INT_MOD_MASK & (generator->r1[hptr] + tmp1);

    new_val = (tmp1 & (~1)) ^ (tmp0 >> 1);

    generator->r0[hptr] = tmp0;
    generator->r1[hptr] = tmp1;

    hptr = --hptr < 0 ? lValue - 1 : hptr;
    hptr = --hptr < 0 ? lValue - 1 : hptr;

    generator->hptr = hptr;
  }

  return new_val >> 1;
}

float get_rn_flt(__global rngen *generator, 
                 int flag)
{
  unsigned int new_val;
  int lptr;
  int hptr = generator->hptr;
  int lval, kval;
  int total_gen = generator->total_gen;
  unsigned int tmp, tmp0, tmp1;
  unsigned int lValue = generator->lValue;
  unsigned int kValue = generator->kValue;

  lptr = (hptr + kValue) % lValue;

  if(flag != 1)
  {
    tmp0 = INT_MOD_MASK & 
      (generator->r0[hptr] + generator->r0[lptr]);
    tmp1 = INT_MOD_MASK & 
      (generator->r1[hptr] + generator->r1[lptr]);

    new_val = (tmp1 & (~1)) ^ (tmp0 >> 1);

    generator->r0[hptr] = tmp0;
    generator->r1[hptr] = tmp1;
    
    hptr = --hptr < 0 ? lValue - 1 : hptr;
    hptr = --hptr < 0 ? lValue - 1 : hptr;

    generator->hptr = hptr;
  }
  else
  {
    tmp = (lptr + kValue) % lValue;
    
    tmp0 = INT_MOD_MASK & 
      (generator->r0[tmp] + generator->r0[lptr]);
    tmp1 = INT_MOD_MASK & 
      (generator->r1[tmp] + generator->r1[lptr]);

    generator->r0[MAP(lptr, total_gen)] = tmp0;
    generator->r1[MAP(lptr, total_gen)] = tmp1;

    tmp0 = INT_MOD_MASK & (generator->r0[hptr] + tmp0);
    tmp1 = INT_MOD_MASK & (generator->r1[hptr] + tmp1);

    new_val = (tmp1 & (~1)) ^ (tmp0 >> 1);

    generator->r0[hptr] = tmp0;
    generator->r1[hptr] = tmp1;

    hptr = --hptr < 0 ? lValue - 1 : hptr;
    hptr = --hptr < 0 ? lValue - 1 : hptr;

    generator->hptr = hptr;
  }

  return (new_val * FLT_MULT);
}

#if 0
int get_rn_dbl(unsigned int gennum, 
               __global rngen *genptr, 
               __global rngen_const *gen_const, 
               int flag)
{
  unsigned int new_val;
  __global unsigned int *r0, *r1;
  int lptr;
  int hptr = gen_const->hptr[get_global_id(0)];
  int lval, kval;
  int total_gen = gen_const->total_gen;
  unsigned int tmp, tmp0, tmp1;

  //lval = gen_const->lval;
  //kval = gen_const->kval;
  r0 = gen_const->bAddr_r0 + gennum;
  r1 = gen_const->bAddr_r1 + gennum;
  lptr = hptr + gen_const->kval;

  if (lptr >= gen_const->lval) 
  {
    lptr -= gen_const->lval;
  }

  if(flag == 1)
  {
    tmp0 = INT_MOD_MASK & (r0[MAP(hptr, total_gen)] + r0[MAP(lptr, total_gen)]);
    tmp1 = INT_MOD_MASK & (r1[MAP(hptr, total_gen)] + r1[MAP(lptr, total_gen)]);

    new_val = (tmp1 & (~1)) ^ (tmp0 >> 1);
    r0[MAP(hptr, total_gen)] = tmp0;
    r1[MAP(hptr, total_gen)] = tmp1;

    hptr --;

    if (hptr < 0)
    {
      hptr = gen_const->lval - 1;
    }

    hptr --;

    gen_const->hptr[get_global_id(0)] = (hptr < 0) ? gen_const->lval - 1 : hptr;
  }
  else if(flag == 0)
  {
    tmp = lptr + gen_const->kval;
    if(tmp >= gen_const->lval)
    {
      tmp -= gen_const->lval;
    }

    tmp0 = INT_MOD_MASK & (r0[MAP(tmp, total_gen)] + r0[MAP(lptr, total_gen)]);
    tmp1 = INT_MOD_MASK & (r1[MAP(tmp, total_gen)] + r1[MAP(lptr, total_gen)]);

    r0[MAP(lptr, total_gen)] = tmp0;
    r1[MAP(lptr, total_gen)] = tmp1;

    tmp0 = INT_MOD_MASK & (r0[MAP(hptr, total_gen)] + tmp0);
    tmp1 = INT_MOD_MASK & (r1[MAP(hptr, total_gen)] + tmp1);

    new_val = (tmp1 & (~1)) ^ (tmp0 >> 1);
    r0[MAP(hptr, total_gen)] = tmp0;
    r1[MAP(hptr, total_gen)] = tmp1;

    hptr --;

    if (hptr < 0)
    {
      hptr = gen_const->lval - 1;
    }

    hptr --;

    gen_const->hptr[get_global_id(0)] = (hptr < 0) ? gen_const->lval - 1 : hptr;
  }

  // new_val = ((unsigned int) temp2 * (double)FLT_MULT + (unsigned int) temp1) * FLT_MULT;

  return (new_val);
}
#endif

/*************************************************************************/
/*************************************************************************/
/*            INITIALIZE: starts the whole thing going                   */
/*************************************************************************/
/*************************************************************************/

CLSprngErr initialize(int ngen, 
                      int gennum, 
                      __global rngen* generator)
{
  int i;
  int j;
  int k;
  unsigned int kValue = generator->kValue;
  unsigned int lValue = generator->lValue;

  //Initialize hptr
  for(i = 0; i < ngen; i++)
  {  
    generator->hptr = lValue - 1;
  }

  si_double(generator->si, 
            generator->nstart, 
            generator->lValue, 
            generator->total_gen);

  get_fill(generator->si, generator->r0, generator);
  generator->si[0]++;
  get_fill(generator->si, generator->r1, generator);

  for(i = ngen - 1; i >= 0; i--) 
  {
    k = 0;

    for(j = 1; j < lValue - 1; j++)
    {
      if(generator->si[j])
      {
        k = 1;
      }
    }   

    if (k == 0) 
    {
      break;
    }

    for (j=0; j < lValue * RUNUP; j++)
    {
      get_rn_int(generator, 0);
    }
  }

  while (i>=0)
  {
    for(j = 0; j < kValue + 1 / 2; j ++)
      get_rn_int(generator, 1);
    for (; j < 4 * lValue; j++)
      get_rn_int(generator, 0);
    i--;
  }   
  return kCLSprngErrNoError;
}

/*************************************************************************/
/*************************************************************************/
/*            INIT_RNG's: user interface to start things off             */
/*************************************************************************/
/*************************************************************************/

CLSprngErr init_rng(int gennum,  
                    int total_gen, 
                    __global rngen *generator, 
                    int seed, 
                    int param)
{
  int doexit = 0;
  int i;
  int length;
  CLSprngErr err;

  if(total_gen <= 0) 
  {
    return kCLSprngErrTotalGenIsNegative;
  }

  if(gennum >= total_gen) 
  {
    return kCLSprngErrGenNumExceedsTotal;
  }

  if(gennum < 0) 
  {
    return kCLSprngErrGenNumIsNegative;
  }

  seed &= 0x7fffffff;  //Remove MSB of seed

  if (param < 0 || param >= kNumParams) 
  {
    return kCLSprngErrUnknownParameter;
  }

  length = validL[param];
  
  //Initialize stuff in the constants part of the generator
  generator->lValue = length;
#if 0
  generator->kValue = validK[param];
  generator->first = validFirst[param];
  generator->total_gen = total_gen;
  generator->LSBS = validLSBS[param];
  generator->seed = seed ^ GS0;
  generator->init_seed = seed;
  
  generator->nstart[0] = gennum;

  for (i = 1; i < length - 1; i++) 
  {
    generator->nstart[i] = 0;
  }

  //err = initialize(1, gennum, generator);
  
  if(err != kCLSprngErrNoError)
  {
    return err;
  }

  
  while(generator->si[0] < total_gen && generator->si[1] == 0)
  { 
    si_double(generator->si, generator->si, length, total_gen);
  }

#endif
  return kCLSprngErrNoError;
}

/*************************************************************************/
/*************************************************************************/
/*                  SPAWN_RNG: spawns new generators                     */
/*************************************************************************/
/*************************************************************************/

int spawn_rng(rngen *genptr, 
              int nspawned, 
              __global int ***newgens, 
              int checkid)
{

  /*
     int **q=NULL, i;
     unsigned *p;
     struct rngen *temp;

     if (nspawned <= 0) 
     {
     nspawned = 1;
     errprint("WARNING","spawn_rng","Nspawned <= 0. Default value of 1 used for nspawned");
     }

     temp = (struct rngen *) genptr;

     p = temp->si;

     q = initialize(temp->rng_type,nspawned,temp->param,temp->seed,p,temp->init_seed);

     if (q == NULL) 
     {
   *newgens = NULL;
   return 0;
   }

   si_double(p,p,temp->lval);

   NGENS += nspawned;

   *newgens = (int **) q;

   if(checkid != 0)
   for(i=0; i<nspawned; i++)
   if(addID((*newgens)[i]) == NULL)
   return i;

   return nspawned;
   */
  return 1;
}

/*************************************************************************/
/*************************************************************************/
/*                  UTILITY ROUTINES                                     */
/*************************************************************************/
/*************************************************************************/
/*
   __device__ int get_llag_rng(void)
   {
   return(lval);
   }

   __device__ int get_klag_rng(void)
   {
   return(kval);
   }
 */

/*
int get_seed_rng(__global rngen *genptr)
{
  return(GS0 ^ genptr->pri.gseed); 
}

int get_hptr_rng(__global rngen *genptr)
{
  return genptr->hptr;
}
*/

//this function return a dynamic buffer create by itself, need further digging to make out how it is used.
__global int* get_fill_rng(__global rngen *genptr)
{
  /*
     int i, *p;
     unsigned int *pp;

     int p[valid[0].L * 2];

     pp = genptr->r0;
     for (i = 0; i < temp->lval; i++)
     {
     p[i] = pp[i];
     }
     pp = genptr->r1;
     for (i = 0; i < genptr->lval; i++)
     {
     p[genptr->lval+i] = pp[i];
     }
   */
  return(NULL);
}

__global int* get_next_index_rng(__global rngen *genptr)
{
  /*
     int i, lval;
     unsigned int *pp;

  //Added by Susan, temprorally use the max length, since dynamic memory allocation is not allowed.
  int p[valid[0].L - 1];

  lval = genptr->lval;

  pp = genptr->si;

  for (i = 0; i < lval - 1; i++) 
  {
  p[i] = pp[i];
  }
   */
  return(NULL);
}

void si_halve(int *a, int length)
{
  int i;

  for (i = 0; i < length - 2; i++) 
  {
    a[i] >>= 1;
    if (a[i+1] & 1)
      a[i] ^= (1 << MAX_BIT_INT);
  }

  a[length - 2] >>= 1;
}

int * get_node_index_rng( int genptr)
{
  /*
     int *p, length;

     length = ( (struct rngen *) genptr)->lval;

     p = get_next_index_rng(genptr);
     if(p == NULL)
     {
     return NULL;
     }

     while (!(p[0]&1)) 
     {
     si_halve(p,length);
     }

     si_halve(p, length);
   */
  return(NULL);
}


/*************************************************************************/
/*************************************************************************/
/*                  MESSAGE PASSING ROUTINES                             */
/*************************************************************************/
/*************************************************************************/


int pack_rng( int *genptr, char **buffer)
{
  return 1;
}


int *unpack_rng( char *p)
{
  return NULL;
}

/*************************************************************************/
/*************************************************************************/
/*      FREE_RNG: remove memory for a generator                          */
/*************************************************************************/
/*************************************************************************/

int free_rng(__global rngen *genptr)
{
  return 1;
}



int print_rng(__global rngen *genptr)
{
  return 1;
}

